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Abstract 

In topology optimization, the design parameter of the element that 
does not give any contribution to the objective function vanishes. This 
causes the stiffness matrix to become singular. To avoid the breakdown 
caused by this singularity, the previous studies employ some additional 
procedures. These additional procedures, however, have some problems. 
On the other hand, convergences of Krylov subspace methods as the so- 
lution of singular systems have been studied recently. By the consecutive 
studies, it has been revealed that conjugate gradient method (CGM) does 
not converge to the local optimal solutions in some singular systems but in 
those satisfying some condition while conjugate residual method (CRM) 
converges in any singular systems. In this article, we will show that a local 
optimal solution of topology optimization is obtained by using the CRM 
and the CGM as a solver of the equilibrium equation even if the stiffness 
matrix becomes singular. Moreover, we prove that the CGM without any 
additional procedures converges to a local optimal solution in that case. 
Computer simulation shows that the CGM gives approximately the same 
solutions obtained by the CRM in the case of a typical cantilever beam 
problem. 



1 Introduction 

In topology optimization problems, we try to obtain the shape and topology of 
the structure that minimizes a given objective function under given constraints. 
The topology of the structure to be obtained in such problems may include the 
number of holes. Although the homogenization method [3] is one of the most 
effective methods for topology optimization in two-dimensional problems, its 
computational cost becomes much larger in three-dimensional problems. For 
this reason, this article treats the material distribution method [2J, which has 
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recently attracted attentions at the computational point of view. In the mate- 
rial distribution method, the shape of the structure is represented as a density 
distribution. The density distribution takes the value one for the element filled 
with material and zero for that containing no material. 

Topology optimization is solved by the following steps ( See FigfT]): 

1. initialize the density distribution of the structure. 

2. repeat the following steps until convergence. 

3. solve the structural analysis problem. 

4. calculate the sensitivity of the density by the sensitivity analysis. 

5. update the density distribution of the structure. 

In this kind of problems, it is always supposed that the design region is divided 
by finite uniform sized mesh, especially the finite element method (FEM) is 
employed in the structural analysis problem. Then the density distribution, the 
displacement field, and the stress field defined on real space are described as 
finite dimensional vectors. 

The densities of some elements of the structure often take the value zero in 
topology optimization. This results in the singular stiffness matrix and the nu- 
merical breakdown of linear solvers in structural analysis. The previous methods 
avoid this breakdown by using the following procedures: 

1. reconstruct the system of equations to regularize the stiffness matrix. 

2. restrict the range of the value of the density pj within < p m in < Pj < 1 
to regularize the stiffness matrix. p mjn is usually set to be 10 -3 in typical 
applications [2J. 

These procedures, however, have some problems. The first one causes the incre- 
ment of the computational cost because the additional steps to find the elements 
with zero density and to reconstruct the system of equations are required. The 
second one causes the physical inconsistency so that the thin or weak material 
exists in the element that contains essentially no material. 

In this article, we will show that such singular systems in topology optimiza- 
tion can be solved by Krylov subspace methods. By these methods without any 
additional procedures, the value of p is not restricted to be zero. Especially we 
are concerned with conjugate residual method (CRM) and conjugate gradient 
method (CGM), the fundamental ones of the Krylov subspace methods. We 
will show that a local optimal solution of topology optimization is obtained by 
using the CRM and the CGM as a solver of the equilibrium equation even if the 
stiffness matrix becomes singular. 

The algorithm using the CGM as a solver of the structural analysis prob- 
lem has already proposed by Fujii et.al. [6] and described its merits for the 
computational cost. However, they used the second additional procedure in ac- 
cordance with the voxel FEM. The convergence of the CGM for the singular 
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stiffness matrix is not considered in their paper. Besides, in CONLIN method 
[4], weli known as a solver of topology optimization problems, the CGM is used 
to obtain the updates of the Lagrange multipliers because singularity 0]. As 
the reason for using the CGM as a solver, they only described that the CGM 
does not need the inverse of the stiffness matrix. No discussion of convergence 
of the CGM to the local optimal solution is also given. 

The convergence of the Krylov subspace methods for a system of linear equa- 
tions with a singular coefficient matrix has been studied recently in numerical 
mathematics. By this consecutive studies, it has been revealed that the CGM 
does not converge to the local optimal solutions in some singular systems but in 
those satisfying some condition [7] and the conjugate residual method (CRM) 
converges to the local optimal solutions in any singular systems [TJ [8] [12] . More- 
over, we will prove that the CGM without any additional procedures converges 
to a local optimal solution in that case. 

The remaining part of this article is organized as follows. First, topology op- 
timization problem is formulated followed by the description of the convergence 
of the Krylov subspace methods. In the section of the Krylov subspace methods, 
we give a brief introduction of the convergence theorems for the CRM and the 
CGM for regular systems and then consider their behaviors for singular systems 
based on the method introduced by Abe et.al. PQ. From this consideration, 
the sufficient condition of the convergence of the CGM for singular systems is 
obtained. After above generic discussion, the convergence of the CGM in topol- 
ogy optimization is proved. Computer simulation for the coat-hanging problem 
verifies our proof. 

2 Formulation of Topology Optimization 

Since topology optimization problems are described in detail in Refs. [2] [5] [6], 
we give the brief introduction here. Besides, we restrict the problem to a two 
dimensional plane strain problem, and suppose to formulate by using FEM for 
comprehensive description. 

The system we are concerned with is, for example, a cantilever with com- 
pletely fixed left end and free right end with a concentrated load at the middle 
point shown in Figj2] The topology and size of the cantilever is determined as 
the solution of the minimization of a given objective function within a given do- 
main. This domain is called design domain and depicted as a shaded rectangle 
in Fig(5] Now we formulate the topology optimization in accordance with the 
Refs. HE]. 

Based on the voxel FEM, a design domain is divided by a uniform sized mesh. 
The material distribution (MD) method is employed for the representation of the 
shape of the structure. In the MD method, the density < p < 1 is assigned to 
each element. The density takes the value zero for the element with no material 
and one for the element with material. A vector composed of all densities is 
called density vector p. Assuming the density directly affects the stiffness of the 
element, Young's modulus is a function of the density. Therefore, every element 
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Figure 1: The flow of computations for topology design using the material 
distribution method and the OC / CONLIN methods for optimization. 
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Figure 2: An example of the topology optimization : a clamped plate. A design 
domain is depicted as a shaded rectangle. 
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of the total stiffness matrix A is a function of the density. The element stiffness 
matrix Aj of the j-th element is also a function of the density pj as follows: 

Ai=f?Pi (1) 

where p is a design parameter often set to be 2 or 3 and Dj is an element 
stiffness matrix seen in ordinary structural analysis problems. This formulation 
is known as Solid Isotropic Material with Penalization (SIMP) model [2J. 

A displacement vector is assigned at every node of the mesh. A rearranged 
vector composed of the displacement vector at every node is called a nodal 
displacement vector afresh and denoted as x. 

The density vector is determined so that the objective function is minimized 
under the given constraints. The objective function is often defined as the mean 
compliance of the structure and the total volume ( or total mass ) and the lower 
and the upper bound of the density are used as the constraints as follows: 

min C(p) := ^ T A(p) X , (2) 
subject to pj < po, (3) 

3 

< Pj < 1 (4) 

The first derivative of the objective function with respect to the design variables 
is called sensitivity. The sensitivity C Pj of Eq.([2]) is given as follows: 

0„ := f. - -A (5) 
dpj d Pj 

Note that x is also the function of p because x is obtained as the solution of 
the linear equation A(p)x = b. Here, b is a vector and the jth element b(j) is 
an external force on the jth node, called nodal force. Considering Eq.(fT]), we 
have 016] 

Pj 

where A^ is the element stiffness matrix of the A:-th element. Xk is defined as 
a vector composed of the elements of the displacement vectors at the nodes 
belonging to the fc-th element. The dimensionality of Xk is, therefore, 8 in 2- 
D space and 24 in 3-D space. The flow of the computation for the topology 
design is shown in FigBfb). The density vector is obtained by the methods for 
optimization problem with constraints, i.e., optimality criteria (OC) method 
and convex linearization (CONLIN) method. Both the OC and the CONLIN 
methods need the sensitivity in Eq.©. The sensitivity is calculated in the 
" Sensitivity analysis" following the calculation of the displacement vector Xj in 
the "Structural analysis" in Fig(5Jb). 

The density of the element that gives no contribution to the objective func- 
tion will vanish. It causes the total stiffness matrix to become singular. In 
such a case, the previous methods avoid the numerical breakdown by using the 
following procedures: 



C Pj = - Xj T (p^)xj (6) 
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1. reconstruct the system of equations to regularize the stiffness matrix. 

2. restrict the range of value of the density to be < p m in < Pj < 1- Pmin is 
usually set to be 10 -3 in typical applications [2J. 

These procedures, however, have some problems. The first one causes the incre- 
ment of the computational cost because the additional steps to find the elements 
with zero density and to reconstruct the system of equations are required. The 
second one causes the physical inconsistency so that the thin or weak material 
exists in the element that contains essentially no material. 



3 Krylov Subspace Methods in Regular Systems 

Before the discussion of the behavior of Krylov subspace method for singular 
systems, we describe the behavior for regular systems. 



3.1 Conjugate Gradient Method 

When a matrix A is positive definite, the solution x* = .<4 _1 b of an equation 
Ax = b is characterized as the minimum point of the following quadratic func- 
tion: 

<t>(x) = ±(x-x*,A(x-x*)) (7) 

= i(x,Ax)-(x,&) + l(x*,Ax*) (8) 

Indeed, because of the positive definitcness of A, <fr(x) > — </>(x*) is satisfied. 
Therefore, by creating the sequence of the vectors XoCxiCx2C- • • which de- 
creases </>(x) , it can be expected to obtain the approximate of the true solution 
x*. This is the essential idea of the CGM. 
The following is one of the CG algorithms: 

CG Algorithm^ 

1. Set the initial guess xo. Similarly, 

r := b-Ax , (9) 
Po := r (10) 

2. For k = 0, 1, • • ■ , repeat the following steps until ||rfe|| < e||b|| 
where e > is a predetermined value: 

( r fc,Pfc) nn 
{p k ,Ap k ) 

Xfc+i = Xfe + QffcPfe, (12) 

r fc+ i = r k -a k Ap k , (13) 

h - -^±i^2, (14) 
(p k ,Ap k ) 

p fc+ i = r k+1 +(3 k p k . (15) 
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Theorem 3-l\g\ 

In the CG method for a symmetry positive definite matrix, the norm 
of the error ||x& — x*|| decreases monotonically. 

3.2 Conjugate Residual Method 

In the CRM, a half of the square of the norm of the residual is used as the 
objective function: 

0(x) = l(Ax-b,Ax-b) = i(r,r) (16) 

and obtain the approximate solution of Ax. = b by successive minimization. 
The following is one of the CR algorithms: 

CR Algorithm^ 

1. Set the initial guess xo. Similarly, 

r := b-Ax , (17) 
Po := r (18) 

2. For k = 0, 1, • • • , repeat the following steps until ||rfc|| < e||b||: 

a k = 7^4^, d9) 
{Apk, Ap k ) 

x k +i = x fc + a k p k , (20) 
r k+1 = r k - a k Ap k , (21) 
(Ar k+1 ,Ap k ) 
(Ap k ,Ap k ) 

Pk+i = r k+1 +(3 k p k . (23) 
The following theorem was proved for evaluating the convergence of the CRM 
Theorem 3-2 

If the symmetric part M = (A + A T )/2 of the coefficient matrix A is 
definite (positive or negative definite ), either of the following holds. 

1. There exists fc s > so that p k ^ 0, (0 < k < fc»-l) and r fctt = 0. 
Further, the following relation holds for < k < — 1: 

K+ill 2 ^ {A mm (M)} 2 (24) 



:|r fe || 2 - X max (A T A) 

is the maximum and minimum eigenvalue of 
the diagonal matrix, respectively. 



x max ^ ''mm 
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2. p k ^ and r fc ^ for all k > and Eq.([2"ip holds. 

The necessary and sufficient condition for the convergence without "break- 
down" of the CRM for an arbitrary b can be derived from the theorem 3-2. 
The derived theorem also provides the meaning of the definiteness of M. Note 
that "breakdown" means "the denominator of the parameter ctk in the conju- 
gate residual algorithm becomes zero so that it becomes impossible to continue 
the computation." The derived theorem said that the CRM converges without 
breakdown for an arbitrary b (for an arbitrary initial guess xo)[T]. 

4 Krylov Subspace Methods in Singular Sys- 
tems 

Now we examine the convergence of the Krylov subspace method for singular 
systems. When dealing with singular systems, the range space R(A) created by 
the row vectors of the coefficient matrix A and the kernel ( or null space ) N(A) 
spanned by x satisfying Ax = play the essential roles. 

According to the analysis by Abe et. al. |T], we prepare the following 
variables: 



where R(A) 1 - is the orthogonal complement of R(A). The following orthogonal 
transformation provides the representation of the Krylov subspace in above 
coordinates: 



where An := QjAQi and A12 := Q\AQ^. By this orthogonal transformation, 
the matrix A can be represented as A in " standard system" where the structure 
of the system is more clear. 

Then, denoting direct sum as ©, the followings hold: 
Theorem 4-1 

An = QjAQx is regular ^ R{A) © N(A) = R n 
Theorem 

A12 = Q1AQ2 = o R(A) 1 - = N(A) 
lemma 4-1 

R(A) 1 - = N(A) An is regular 



r := rank(A) = dim(R(A)) > 0, 

qi , ■ • • , q r : orthonormal basis of -R(^4) 

q r +i, • • • , q« : orthonormal basis of R{A)- 

Qi '■= (Qi) ' ' ' > <3r) ■ n x r matrix 

Q2 ■= (qr+i) 1 ' • j Qn) :nx(n — r) matrix 

Q := (Qlt Q2) '■ n x n orthogonal matrix 




(25) 



9 



From theorem 4-2, when the condition " R(A) 1 - = N(A)" is satisfied, the 
standard system reduces more simple structure as follows: 

A = Q T AQ=^ °) (26) 

and, from lemma 4-1, An = Q\AQ\ is regular in this case. In topology opti- 
mization, the partial stiffness matrix is regular even if the total stiffness matrix 
is singular, we can use Eq. (|2"6")l as the standard system. 



4.1 Convergence Property of CRM 

It is guaranteed for the CRM to converge to the local optimal solution for a 
system of linear equations with a singular coefficient matrix jl D For comparison 
with the CGM described later, we show the method of analysis described in [T]. 
The coefficients a k and (3k in the CRM is calculated as follows: 

(r fc , Apk ) 

ak - 



fik = 



(Ap k , Ap k )' 
{Ar k+1 ,Ap k ) 
(Ap k , Ap k ) 



Using the standard system Eq.(|26|). these coefficients are decomposed into the 
R(A) component (with superfix ||) and N(A) component (with superfix _L), 



A n \ pi 

)\ pi 



An V PM M" U P 

/ 

(4^hp!) 



(i l! M pi )'{ o o ) { pi 



(A u pl,A u pl) 



An \ / P H An \ p" 

vUJ'l J\ pi 

(^ur| +1 ,A llP |) 
(Anpl,Anpl) 

Using above equations, the algorithm of the CR method can be decomposed 
into the R(A) and the N(A) component. 



The decomposed CR algorithm 
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1. Choose x{J and Xq. Similarly, 

R(A) component N(A) component 

r D -^-11 X 0! r o — ° > 

n ll _ -II n -L _ _-L 

Po — r o> Po — r 

For 

2. k = 0, 1, ■ • ■ , repeat the following steps until the R(A) component of the 
residual converges: 



R(A) component -^(^4) component 



(AupJUupJiy 



x fe+i = x fe +"fcPfc, x fc+i = x fe + afcPfe , 

r l+i = r|i + "fc^nPL r fc+i= r fc> 
_(r| ±11 Anp|) 

^ (pE^up!) ' 

p!+i = r I+i + &p|. Pfe+i = r fe+i + /^fcPfe 



Since the algorithm for the R(A) component can be regarded as the CRM 
applied to the R(A) subsystem, it is guaranteed that the norm of the R(A) 
component of the residual decreases monotonically from theorem 3-3. Since the 
N(A) component of the residual is equal to b- 1 and unchanged, consequently for 
an arbitrary b the convergence of the CRM for singular systems is guaranteed. 
On the other hand, there has no information about x -1 . However, when h 1 ^ = 0, 
r 1 - = and p- 1 = and then x^ is constant. 

Because of the convergence of the CRM in singular systems as described 
above, we can use the CRM to obtain a local optimal solution in topology 
optimization even if the stiffness matrix becomes singular. In terms of compu- 
tational cost, however, the CGM will be much effective than the CRM. Then 
we will examine the convergence of the CGM in singular systems. 

4.2 Convergence Property of CGM 

Similar to the previous section, we look the CGM in the standard system below. 
The fundamental steps of the CGM are described in the previous section. 
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For analyzing the behavior of the R(A) component and the N(A) component 
of the CGM, first we decompose the vectors x, p, b, and r into these subspaces. 



X 


= Q T x = 


(xll 




p 


- Q t p = 


(P 1 


l ,F 1 ) T 


b 


= Q T b = 


(bl 




f 


= Q T r = 


(r". 





Eqs. (|12ll3ll5[) can be decomposed into the R(A) and the N(A) component as 
follows: 

4+1 = x li + Qfc Pfe' x fc+i = x ^ + "^Pfc ' 

pH+i = Pfc + Abp|, Pfc+i = Pfc + QfcPfe • 
Then otk and /3& can be rewritten as follows: 

.0 \ / p | 
Pi 



p! W w p^ 



/ v j \ pt 

.ii 



(pjl.AnpJi) 



/3 fe = 



r ll 

_JL 
r fc+l 



pI 

Pi 







(pjUnpjl) 




' A o !)( 


'p! 




v Pfc 


An 0\/ 


pJ! 


o o ; 


Pfc 1 



(27) 



(rl +1 ,A llP l) 
(pL^iipI) 



(28) 



Using above descriptions, we try to decompose the CG method into the R(A) 
and the N(A) component. 



The behavior of the R(A) and the N(A) component of the CGM 

1. Choose xjj and Xq-. Similarly, 

R(A) component N(A) component 



b" -A n x|l, 
p8=rl, 





= b x 




- r 1 - 


Po 


— r o 
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2. For k = 0, 1, • • • , repeat the following steps until the R(A) component of 
the residual converges: 

R(A) component N(A) component 



o 



II _ rf.pJl) ■ (4, Pi) 



k (pI.AupI)' " (pf^iipl)' 

4+i = 4 + ^pL x fc+i = x fe +«fePfe . 

l+i = r\ + a fc ^iip|, r£ +1 = r£, 

a ( r fe+l = A HP|) 

Pfe = 



Pfe+i = 4+i + hvl Pfe+i = ri +1 + f3 k p£ 



For the CGM for symmetric definite matrix, it is guaranteed that the norm 
of the error decreases monotonically from theorem 3-1. Therefore, if the R{A) 
component of the CG algorithm is closed in the R{A) subspace, it is guaranteed 
that the norm of the error decreases monotonically. Since the N(A) component 
of the residual is equal to b and unchanged, the convergence of the CGM 
for singular systems is guaranteed in such a case. However, focusing on the 
coefficient a^, as described in Eq. ([271 . we can see its numerator includes 
and Pj7, the N(A) components of rfc and p k - Since the R(A) component of 
the CGM is not closed in the R(A) subspace, the convergence of the CGM for 
singular systems is not guaranteed. Indeed, the numerical experiments that the 
CGM diverges when b^ ^ has been reported [7]D If b G R(A), that is h 1 - = 0, 
then rjf — 0, — 0, (k — 1, 2, • • • ), and aj: = 0, (k = 1, 2, • • • ), resulting in the 
convergence of the R(A) subspace. This is the sufficient condition of convergence 
of the CGM for singular systems. 



5 Singular Systems in Topology Optimization 

When the densities of all elements adjacent to a node j takes the value 0, all 
the corresponding components of the total stiffness matrix vanish. Then the 
total stiffness matrix becomes singular and numerical algorithms for solving the 
equilibrium problem might be break down. In order to make the matrix regular, 
reconstructing the system of linear equations by removing the j-th row and the 
j-th column. Since the reconstructed stiffness matrix is considered to be An 
in Eq. ([26| , the range space of the reconstructed stiffness matrix corresponds to 
R(A) and the removed part N(A). If the j-th component of the nodal force 
vector b takes the value 0, then b e R(A), and from the discussion in the 
previous section it is guaranteed that the CGM converges without breakdown. 
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Now the proposition to be proved as the sufficient condition for the conver- 
gence of the CGM is as follows: 
Proposition- 1 

If the densities of all elements adjacent to the node j takes the value 
0, the nodal force at the node j is 0. 

To prove the proposition- 1 directly is difficult because the proposition- 1 must 
evaluate the nodal forces after giving the densities of the elements while the 
practical computational process calculates the densities after giving the nodal 
forces. Then we prove the contraposition of the proposition- 1: 
proposition-2 

If the nodal force at the j-th node is not 0, there exists at least one 
element which density takes the positive value among those adjacent 
to the element j. 



5.1 Proof of Proposition-2 

Assuming b(j) ^ for Ax = b, then 

k—rij 

? b(j) = :)x = J2 Mik, :)a* (29) 
fc=l 

where k is an index of an element adjacent to a node j, and Ak are the 
element displacement vector and the element stiffness matrix of the k-th element, 
respectively, jj. is the local index of the j-th node among the nodes located on 
the boundary of the fc-th element. Ak(jk, '■) is the row vector corresponding to 
the jfc-th node of Ak- Among nj terms in the left hand side of above equation, 
there exists at least one element, say the k-th element, so that x-^ ^ and 

%(4> : )x ¥ ^o. 

The sensitivity Cp_ of C(p) with respect to the density of the fc-th element 
is given from the eq.©: 



^ = -4(p^K ( 3 °) 



For k, since A^ is symmetry positive definite, Cp^ < 0. 

Next, we evaluate the absolute value of C^. Using x = A~ x b and Eq.([T]), 



Cp- can be rewritten as follows: 



T °k U k "k 

Since p, thr, and Dx are constant, we have 



= -pfC^D^ (31) 



lim Co- = -oo (32) 
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Based on the above, we will show that the density px of the fc-th element 
adjacent to the j-th element ( bj ^ ) moves away from the value even if 
the initial guess of px is close to 0. px is obtained as a solution of constrained 
optimization problem by iterative methods in general. Then we examine for the 
following three principal methods: 

1. gradient vector based method including steepest descent method and CGM 

2. the OC method 

3. the CONLIN method 

First, for the gradient vector based method, since px is updated to the 
negative direction of Cp^, px definitely increases and does not converge to zero. 

Next, we examine for the OC method. The updating equation of px is as 
followsgjF 



, ( t ) v 0.85 



A 

where A < is a Lagrange multiplier. From the consideration of the case when 
px is close to zero, there exists px > so that — A > 1. Since at that point 

Px > Px\ the above updating equation moves px away from 0. 

The updating equation in this case by the CONLIN method is derived from 
ref.[4] as follows: 

4' +1) = (-^) 1/2 4" <«> 

where A > is a Lagrange multiplier. Similar to the OC method, From the 
consideration of the case when px is close to zero, there exists px > so that 

( > 1. Since at that point p^ +1 ^ > px \ the above updating equation 

moves px away from 0. 

As described above, px does not converge to zero when solving by the iter- 
ative methods with the initial value of < px < 1- Now the proposition-2 was 
proved and so was proposition- 1. But it must be noted that if pk is initialized 
to be 0, the CGM breaks down at the time when calculating its sensitivity C Pk . 



6 Simulation 

What is verified in computer simulations is that the CGM converges without 
breakdown to a local optimal solution even if the total stiffness matrix becomes 
singular. 

One can see, however, that the convergence rate of the density to the value 
is much slow while the rate to the value 1 is very fast both in the OC and the 
CONLIN method. That is, the density going toward zero takes the tiny values 
such as 10 _ln , 10~ 20 , and 10~ 30 endlessly. Then the total stiffness matrix 
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does not easily become singular. Therefore, we add the operation to force the 
density value be if it is less than a predetermined value ( 1CP 3 in the following 
simulations ). Besides, we use a simple pre-conditioning matrix for CG and 
CR method, every diagonal element of which are equal to the inverse of the 
corresponding diagonal element of the stiffness matrix. Pre-conditioned CG 
and CR is denoted as PCG and PCR, respectively. 

We compare the solutions obtained by the following four methods: 

1. PCG-OC 

2. PCR-OC 

3. PCG-CONLIN 

4. PCR-CONLIN 

The environment of the computer simulation is as follows: 

• CPU: PentiumV(933MHz) 

• memories : 512MB 

• compiler : gcc version 2.95.3 20010315 (release) 

• compiler option : 02 

Some pre-determined values are as follows: 

• upper bound of the total volume of the structure: pmax / X^" Pi ~ 0-375 

• convergence criteria of structural analysis problem: repeat until either the 
residual becomes less than 10 -8 , or the number of iterations exceeds a 
predetermined value. The predetermined value is set to be the number of 
nodes. 

• convergence criteria of the OC or CONLIN method: repeat until either the 
absolute value of the change of Lagrangian becomes less than 10 -10 , or the 
number of iterations exceeds the predetermined value. The predetermined 
value is set to be 100. 

6.1 Two bar truss problem 

A simple problem for which the solution is obtained analytically is considered as 
a verification of the previously described proof. Two bar truss problem is well 
known one of such problems. Figl^a) illustrates the problem definition. The 
design domain, a — 10 [mm], b — 20 [mm], is discretized using a 20 x 40 mesh of 
four-node bilinear plane strain elements. The material parameters are assumed 
to be Young's modulus E = 2.1 x 10 5 [7V/mm 2 ] and Poisson's ratio v — 0.3. The 
load at the middle of the free end is assumed to be P = 1.05 x 10 2 [iV]. 
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( b ) Final topology layout obtained ( c } change of total strain energy 

by PCG-CONLIN and PCR-CONLIN. in the process of optimal solution search 



Figure 3: simulation results for two bar truss problem 
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method 


total iterations 


total strain energy 


CPU time (sec) 


PCG-OC 


2599 


0.0191866 


141.66 


PCG-CONLIN 


2014 


0.0160593 


111.78 


PCR-OC 


45031 


0.0191104 


1958.74 


PCR-CONLIN 


40033 


0.0160725 


1590.87 



Figure 4: Comparison between methods for two bar truss problem 



The final topology layouts using the four methods are given in FigEfa) and 
(b). All solutions give so-called two bar trusses with an internal angle of 90° 
which is exactly the same as the analytical solution [TQl [11] . 

FigEJc) shows the total strain energy iterative histories for the four meth- 
ods. The horizontal and vertical axis in Fig|3tc) is the cumulative number of 
iterations of PCG or PCR methods and total strain energy, respectively, in log- 
log scale. In the histories, we can see that PCG-CONLIN is the fastest among 
others. Indeed, as denoted in Fig|4] PCG-CONLIN obtains the minimum total 
strain energy by the minimum CPU time among others. In the previous part of 
this article, we showed the convergence of the CRM and the CGM for singular 
systems. The computer simulation verified that both of them can use for the 
structural analysis in topology optimization even if the stiffness matrix becomes 
singular. It also showed that the CGM is much effective than the CRM in terms 
of computational cost. 

7 Conclusion 

In topology optimization, a singular stiffness matrix is often encountered be- 
cause the densities of some elements of the structure become zero. To avoid 
the numerical breakdown caused by this singularity, the previous methods use 
some additional procedures for regularizing the singular matrix. For example, 
the lower bound of the density was introduced for the material distribution 
method. These procedures have, however, some problems. To resolve such a 
singular system without additional procedures, we focused on the convergence 
properties of the conjugate residual method (CRM) and the conjugate gradient 
method (CGM). The convergence of the CRM for singular systems has been 
proved but has not been applied to the structural analysis in the previous work. 
In this article, computer simulation showed that using the CRM as a solver of 
the structural analysis gives a local optimal solution in topology optimization. 
Next, the convergence of the CGM for singular systems, especially in the case 
when the stiffness matrix becomes singular in topology optimization, considered 
in this article. Although the idea that uses the CGM for solving the structural 
analysis in topology optimization has already been proposed, the proof of its 
convergence has not been given. We proved the convergence of CGM when the 
stiffness matrix becomes singular. It also holds for preconditioned CGM. Com- 
puter simulations for an analytically solved problem verified our proof. Because 
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of no restriction for design variables, the similar discussion can be available for 
the homogcnization method. 
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